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Abstract 

We consider the linear stability of shear banded planar Couette flow of the Johnson- 
Segalman fluid, with and without the addition of stress diffusion to regularise 
the equations. In particular, we investigate the linear stability of an initially one- 
dimensional "base" flow, with a flat interface between the bands, to two-dimensional 
perturbations representing undulations along the interface. We demonstrate ana- 
lytically that, for the linear stability problem, the limit in which diffusion tends 
to zero is mathematically equivalent to a pure (non-diffusive) Johnson-Segalman 
model with a material interface between the shear bands, provided the wavelength 
of perturbations being considered is long relative to the (short) diffusion lengthscale. 

For no diffusion, we find that the flow is unstable to long waves for almost all 
arrangements of the two shear bands. In particular, for any set of fluid parameters 
and shear stress there is some arrangement of shear bands that shows this instability. 
Typically the stable arrangements of bands are those in which one of the two bands 
is very thin. Weak diffusion provides a small stabilising effect, rendering extremely 
long waves marginally stable. However, the basic long-wave instability mechanism 
is not affected by this, and where there would be instability as wavenumber k — > 
in the absence of diffusion, we observe instability for moderate to long waves even 
with diffusion. 

This paper is the first full analytical investigation into an instability first docu- 
mented in the numerical study of [1]. Authors prior to that work have either hap- 
pened to choose parameters where long waves are stable or used slightly different 
constitutive equations and Poiseuille flow, for which the parameters for instability 
appear to be much more restricted. 

We identify two driving terms that can cause instability: one, a jump in N±, as 
reported previously by Hinch et al. [2]; and the second, a discontinuity in shear 
rate. The mechanism for instability from the second of these is not thoroughly 
understood. 
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We discuss the relevance of this work to recent experimental observations of 
complex dynamics seen in shear-banded flows. 

Key words: shear-banding fluid, linear instability, Couette flow, diffusive 
Johnson-Segalman fluid, interfacial instability 
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1 Introduction 



Complex fluids such as wormlike [3] and onion [4] surfactants commonly show 
flow instabilities and flow-induced transitions that lead to spatially hetero- 
geneous, "shear banded" states. In shear thinning wormlike micelles, for ex- 
ample, homogeneous flow becomes unstable above a critical shear rate. The 
system then separates into bands of differing viscosity and internal structure, 
separated by an interface that has its normal in the flow-gradient direction. 
Widespread experimental observations of this phenomenon have been made 
by flow birefringence [5], and by NMR [6] and ultrasound velocity imaging [7]. 
More recently, fluidity banding has been reported in soft glassy materials [8]. 
In bulk mechanical measurements, the main signature of shear banding is a 
kink followed by a plateau in the steady state flow curve. 

Beyond this basic picture, an accumulating body of data reveals that shear 
banded states can fluctuate. Early evidence came from unsteady erratic [9] 
or periodic [10,11] fluctuations in the wall stress at an applied value of the 
shear rate. More recent velocimetry experiments with enhanced spatial and 
temporal resolution have unambiguously revealed fluctuations in the interface 
between the bands [7,12,13,14,15]. To date, however, most theoretical studies 
have considered only a flat, stationary interface. In this paper, therefore, we 
study analytically the linear instability of shear banded flow with respect to 
small undulations along the interface. 

Theoretically, shear banding is thought to arise from a non-monotonicity in 
the underlying constitutive relation between the shear stress and shear rate 
for homogeneous flow [16,17,6,18]. The simplest constitutive model to mimic 
this dependence (apart from "toy" models that do not obey the principle 
of material frame indifference) is the Johnson-Segalman (JS) model [19]. A 
sample plot of shear stress against shear rate in simple shear flow for this fluid 
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Fig. 1. A typical plot of shear stress against shear rate for steady, homogeneous 
shear flow of a shear-banding fluid. Here we use the Johnson-Segalman model with 
parameters e = 0.05, a = 0.8. If the shear stress is T&, homogeneous shear flow with 
any of the shear rates jl < 7m < in is permitted; but flow with shear rate 7m is 
unstable to one-dimensional perturbations. 



is given in figure 1. Homogeneous flow with a shear rate on the decreasing 
part of the curve (e.g. 7m) is unstable to one-dimensional perturbations with 
wavevector in the flow-gradient direction [20]. The fluid therefore separates 
into a structure comprising bands of differing shear rates 7l and jn, one on 
each of the stable, upward-sloping parts of the curve. The interface between 
the bands has its normal in the flow-gradient direction. The shear stress T is 
uniform across the whole flow, as required by a force balance. 

The JS model in its original form contains no mechanism for uniquely selecting 
the shear stress T b at which banding occurs. Instead, in a numerical study of 
shear banding in such a model, the shear stress in the steady banded state 
depends strongly on the startup history [21,22,23,24], and can have a stress 
anywhere in the range Ti < T b < T 2 in figure 1. This conflicts notably with 
experiment, which consistently reveals a highly reproducible banding stress. 

It is therefore critical to regularise the model in some way, to ensure stress 
selection. This is achieved by modifying the constitutive equation to include a 
diffusive ("non-local") term [22,25,26,27]. This mechanism was first proposed 
in 1989 by El-Kareh & Leal [28], with the physical interpretation that individ- 
ual polymer molecules can slowly diffuse across the interface, carrying their 
stress histories with them. Such terms also arise naturally in models of liquid 
crystalline dynamics. Regardless of their physical origin, non local terms lead 
to the selection of a unique, reproducible shear banding stress, T b , as seen 
experimentally. They also provide a length scale for the thickness of the in- 
terface. In contrast, in the local model the interface is unphysically sharp: the 
flow variables jump discontinuously across it. 

The exact details of how stress diffusion should be added to the JS model vary 
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from author to author. The most commonly used version is due to Olmsted & 
coworkers [29,30], in which a term proportional to the Laplacian of the extra 
stress, V 2 S, is added to the evolution equation for the polymeric stress. In 
this paper we will use a slightly more general form that is also capable of 
incorporating the model introduced by Yuan [31], in which the term added is 
a negative multiple of the Laplacian of the rate of strain. 

The theoretical framework just described has been developed largely in the 
context of one dimensional (ID) studies that consider only the flow gradient 
direction, normal to the interface between the bands [21,32,33]. Clearly, such 
studies assume from the outset that the interface between the bands is per- 
fectly flat and they predict (with few exceptions: [34,35]) time- independent 
banded states. This is clearly at odds with the accumulating body of data 
described above, revealing fluctuations of the banded state. 

In view of this, a crucial question is whether the stationary, flat banded state 
of ID calculations will persist in 2D, or whether it destabilises to exhibit large- 
amplitude interfacial fluctuations. The first step to answering this is clearly to 
perform a linear stability analysis of the ID "base state" with respect to small 
2D (flow, flow-gradient) perturbations corresponding to wavelike undulations 
along the interface. 

This was first addressed within local models. McLeish [36] considered a Doi- 
Edwards type fluid in capillary flow. He found instability to long waves, pro- 
vided the high shear rate band is very narrow. As we shall see later, this is 
qualitatively very different from our results. (McLeish did not give the spe- 
cific parameters of his calculation, so quantitative comparison is not possible.) 
Renardy [37] examined the stability of the local JS model in planar banded 
Couette flow. She found linear instability for short wavelengths (wavenumber 
greater than 8). For mainly historical reasons, however, she happened to con- 
fine her study to a base state corresponding to "top-jumping" (T b = T 2 ) and 
an extremely thin high-shear band. We will return below to comment on this 
choice in the context of our own findings. 

The first observation of linear instability within the diffusive JS model was 
in the numerical study of [1]. This considered general band thicknesses and 
demonstrated, for the first time, instability with respect to long and moderate- 
length waves. The short-wave instability predicted in [37] was eliminated by 
the stabilising presence of diffusion. A subsequent non-linear numerical study 
showed the interface to be restabilised at the level of finite amplitude fluctu- 
ations [38]. 

The main contribution of the present study is a detailed analytical interpre- 
tation of the numerical findings of [1]. We start by deriving the important 
result that, for the limit of weak diffusion, the 2D linear stability properties 
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of the diffusive JS model are equivalent to those of the original local model, 
with a "material interface" (defined below) between the bands. This equiva- 
lence was not obvious a priori, since the ID behaviour differs so dramatically 
between the two models: as noted above, the diffusionless limit is singular in 
ID because it has no mechanism for selecting the banding stress. The addition 
of weak stress diffusion thus drastically modifies the ID global properties, by 
selecting a unique base state out of the continuum of possibilities. 

This equivalence with the local model allows us to simplify considerably the 2D 
stability analysis of the diffusive case. In consequence, we are able to plot out 
the full spinodal boundary of instability in the phase diagram, and to predict 
the dispersion relations seen numerically. For a very wide range of model 
parameters, we find instability to waves of moderate wavelength A having 
h 2 <C A 2 < L 3 /h for small diffusion length, h (where L is a typical channel 
lengthscale) . We also identify two possible driving terms for instability: one 
due to a jump across the interface in the base state shear rate, the other due 
to a jump in base state normal stress. 

The paper is structured as follows. In section 2 we introduce our governing 
equations; in §3 we lay out what is known about the one-dimensional steady 
state solution of these equations. In §4 we set up the two-dimensional linear 
stability problem, and in §5 show analytically that, in the limit of small dif- 
fusion, for perturbations whose wavelength is not asymptotically short, the 
region between the shear bands may be considered as a material interface. 
This allows us to work, for the remainder of the paper, with the simpler non- 
diffusive model (but with a value for the selected stress in the one-dimensional 
base state, as selected by diffusive terms). A long- wave stability analysis 
is given in §6, and the results of a full numerical calculation in §7. In §8 we 
review our results and draw conclusions. 



2 Governing equations and dimensionless form 

The standard equations governing flow of an incompressible inertialess fluid 
are conservation of mass and a force balance: 



in which U_ is the velocity field and S_ is the total stress tensor. The stress 
consists of an isotropic pressure P, a Newtonian solvent term of viscosity rj 
plus a polymer extra stress 



V • U = 



V -S = 




S = - 



PL + ^r\K + 



(2) 
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The extra stress evolves with dynamics based on the Johnson-Segalman equa- 
tion but with two possible diffusion terms added to regularise the equations: 



A I 

2GR - -g - 2TV 2 K + VV% (3) 

The case T> = reproduces Yuan's model [31]; the case T = reproduces the 
more usual model of Olmsted & coworkers [29,1]; and if V = T = we regain 
the original Johnson-Segalman model [19]. 

The derivative is the Johnson-Segalman form, with "slip parameter" a: 
A ( 9 \ 

^ij— I t^- + U k V k J — [H ik Q k j — Qik^kj + a{Ei k Y, k j + T, ik E k j)} , (4) 
and the tensors E_ and O are based on the velocity gradient: 

Eij = \{ViU 3 + VjUi) n tJ = |(V,^ - V^)- (5) 

We scale lengths with the channel width L, times with the polymer relax- 
ation time r, and stresses with the modulus G. The dimensionless governing 
equations are then: 

V-C/ = V •£ = (). (6) 
g = -PI + 2eR + ^ (7) 

A , . 

g= 2^-^ + / 2 {-2AV 2 ^+ (1 - A)V 2 ^}. (8) 
We have used three dimensionless parameters: 

e = v /Gr l 2 A = TG/L 2 l 2 (l - A) = Vr/L 2 . (9) 

Thus e is the retardation parameter, / is a characteristic diffusion lengthscale, 
and A is a selection parameter to choose between the different stress-diffusion 
mechanisms. The dimensional diffusion length h is IL. Shear-banding is pos- 
sible for e < 0.125 and in all our examples we will use the value 0.05; it is 
more difficult to predict the likely physical values of other parameters, so we 
consider ranges < I < 0.01, 0.3 < a < 0.8 and < A < 1. 

We will be considering a planar shear flow, bounded between two solid walls. 
As well as the standard velocity boundary conditions of no slip and no pene- 
tration at the walls, we require boundary conditions on the stress because of 
the higher derivatives introduced by the diffusion terms. We have chosen to 
use a boundary condition imposing no flux of extra stress at the boundaries: 
for a wall given by y = constant this imposes the condition <9S/ dy = at the 
wall. Physically this constrains our steady flow not to have boundary layers 
near the walls, although (as in §3.2) boundary layer structures are permitted 
within the body of the fluid. 
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3 Steady-state solution 



If we impose the constraint of a steady unidirectional flow with no pressure 
gradient along the channel, then all flow variables depend only on the position 
across the channel, y. We denote differentiation with respect to y by D. We 
also introduce new variables T^- based on the stress tensor E^-: 



T 11 = i(l-o)E 11 -i(l + a)E 22 (iO) 
T 22 = ±(l-a)E n + ±(l + a)E 22 (If) 
Ti 2 = S 12 . (12) 

The equations governing the flow then reduce to 

U x = U(y) U y = 7 = DU (13) 

for the velocity field, and 

S n = -P + (f - a)-\T 22 + T n ) (14) 
Sn = ej + T 12 (15) 
S 22 = -P+(l + a)-\T 22 -T 11 ) (16) 

for the total stress. Finally, for the polymer stress we have 

Z 2 (1-A)/J 2 T 22 -T 22 = (17) 

which is satisfied if T 22 = 0, and 

(f - a 2 ) 7 T 12 = -Z 2 (f - A)D 2 T n + T n (18) 
7T11 =7 - T 12 - / 2 AD 2 7 + Z 2 (l - A)D 2 T 12 . (19) 

In Couette flow, there is no pressure gradient along the channel: 

= 0, (20) 

and the unsealed momentum equations are 

T 12 + e 7 = T b (21) 

where we recall that the shear stress T b is selected from the continuum of 
possibilities Ti < T b < T 2 only in the diffusive model [29] , and 

dP 

= -(l + a )~ 1 DT ll . (22) 

dy 
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This last equation is the only place that the dependence of P on y appears, 
so we can solve it by setting 



Now let us look at two cases: / small but non-zero (diffusive JS), and I = 
(non- diffusive). If / = we have the Johnson-Segalman model in its origi- 
nal form, with no mechanism to select the shear stress: in this (unphysical) 
case, shear banding can occur at any shear stress in the interval 7\ < T b < 
T 2 (see figure 1). For any such value, there are three possible shear rates 
7l(2&) < 7m (Tb) < 7h(7&), of which a flow with would be unstable to 
one-dimensional perturbations. The system will therefore consist of bands of 
the two shear rates tl and 7^. There may be several of these. 

If we allow / 7^ 0, the system changes radically. The shear stress T5 = T sc \ 
is uniquely selected [39,29], and the flow will separate into shear bands of 7^ 
and of 7/f. Between these bands, rather than a sharp interface, is a matching 
region of width /, across which the flow variables vary continuously. 

In order to carry out direct comparisons between the / = and I 7^ cases, we 
will assume here that only one region of each shear rate is formed, and that 
there is then only one interface between high- and low-shear rate bands. 

For definiteness, we assume that there is a band of low shear rate 7^ near 
the wall y — 0, which continues up to a position y — k. Near the wall y = 1 
(and continuing down towards y — k) there is a band of high shear rate jn- 
The average shear rate (and hence the speed of the upper wall) is U wa \\ = 
+ (1 — K )lH- This can equally be considered as the Weissenberg number in 
our nondimensionalisation. The symmetry of planar Couette flow is such that 
these two bands could be interchanged to produce an essentially equivalent 
flow: we restrict our attention to the arrangement with the low-shear band 
near y — 0. 

For 1^0, the two regions well away from the matching region will be denoted 
the outer, and the matching region of width /, the inner region. Each outer 
region is undergoing homogeneous shear flow, while the inner region can be 
thought of as a finite-width interface between these two phases. 

3. 1 Outer solution 

Well away from the matching region, we expect derivatives of all our quantities 
to be at most order 1: so if I is very small we can neglect terms l 2 D 2 and 
regain the equations governing one-dimensional flow of the original Johnson- 




(23) 
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Segalman equation. This system has been investigated thoroughly [29,33]. The 
stress components of the solution are 

(l-a 2 )7 2 7 
Tl1 = 1 + (1 - a^ 2 Tl2 = 1 + (1 - a^ 2 T22 = ° (24) 

and the shear stress condition T± 2 = — gives a cubic equation for 7: 

(1 _ a 2 )e ^3 _ (1 _ a 2 )TbAf 2 + (1 + £ ^ _ Tfc = q (25) 

Since the neglected terms all involve derivatives of the leading-order terms, 
which are constants, this is an exact solution to the governing equations for 
each permissible value of 7. 

If a 2 < 1 and e < 1/8, there is a range of values of the shear stress TJ, for 
which Ti < Tf, < T 2 and equation (25) has three solutions tl, 7m and jn, 
which allows the possibility of shear bands. 



3.2 Inner region 



If I is small but non-zero, we do not expect a truly sharp interface between 
regions of high and low shear rates. Instead, in the inner region, the second- 
order derivatives become 0{l~~ 2 ) and so become large enough to be important, 
despite being multiplied by I 2 . 

We take this region to be a layer centred on k and use a rescaling variable 
y = k + /£. The equations become (denoting derivative with respect to £ by 
d): 



(1 - a 2 ) 7 T 12 = -(1 - A)rf 2 T n + T n (26) 
7 T n = 7 - T12 - Ad 2 7 + (1 - A)d 2 T 12 (27) 
T 12 + e 7 = T 6 (28) 

with matching conditions at the edges of the layer for 7, T n and T 12 in terms of 
the outer solution. Existence of a solution to this system which matches onto 
the two constant-7 solution branches (with 7 dependent on T b ) as £ — > ±00 is 
a constraint on T b . 

This is most easily seen in the case A = 1, in which, following Lu, Olmsted & 
Ball [25], (28) can be written as 

rf 2 7 = 7 + ( e7 -T 6 )(l + (l-a 2 ) 7 2 ). (29) 
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Multiplying by d'j and integrating with respect to £ across the interface region 
gives a continuity condition 



](1 + e)f + " « 2 )7 4 - T b (7 + ±(1 - a 2 )7 3 ) 

which, coupled with the fact that 7^ and 7# depend on Tb, has a unique 
physically relevant solution 

T b (e,a) = ^(l-a 2 r^(2e-Ae 2 ) 1/2 - (31) 

A similar result (with a different function TJ,(e, a)) is true for other values of 
A. For reference, we note that at e = 0.05, the extremal selected values are 

A = 0: T h = 0.48284(1 -a 2 )" 1/2 ; A = l: T b = 0.45(1 - a 2 )~ 1/2 . (32) 

In this way, the addition of weak diffusion is seen to regularise the equations 
by selecting one possible value T b out of the continuum of possibilities 7\ < 
Tb < T2 at which a one-dimensional shear-banded solution can occur. This 
one-dimensional solution will be taken as the relevant base state about which 
we analyse two-dimensional perturbations (in the flow, flow gradient plane) in 
the remainder of the paper. 

4 Linearised equations 

We now add a small perturbation proportional to £ , such that 

LL=(u + uS,vs) P = P hasc +p£, (33) 



= (30) 



S 



S\l + S\\£ S12 + s^S 



S\2 + s 12^ S22 + S22& j 



T 



Tu + tu£ T12 + t 12 £ 

+ t±2£ ^22^ 



(34) 



in which 

£ = <5exp (ikx — iut) (35) 

for a small parameter 5, and we will use D to denote differentiation with 
respect to y. We have introduced the dimensionless wavenumber k = 2nL/\ 
for waves of wavelength A. We can now linearise the equations about the 
base state in this small parameter, to obtain our new governing equations for 
evolution of the linear system. For a fixed real wavenumber k, if we obtain 
a valid solution with lm(u) > then the system is linearly unstable. This 
analysis is necessarily restricted to considering instabilities which manifest 
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at linear order; any instability (like the high-Reynolds number instability of 
laminar flows) which appears only at nonlinear order will not be captured 
here. The mass conservation equation is 



iku + Dv = 0. (36) 

This is automatically satisfied by the introduction of a streamfunction ip with 
v = —ikip and u = Dip, such that the perturbation velocity uS = V x (ipz£). 
The remaining governing equations for the force balance become 

iksu + Dsi2 = iksw + Ds 22 = (37) 

sn = -p + 2dkDiP + (1 - a)-\t 22 + t u ) (38) 
s 12 = e(D 2 4> + k 2 iP) +t 12 (39) 
s 22 = -p- 2eikDiP + (1 + a)-\t 22 - t n ) (40) 



and for the polymer stress we have: 

{-iuj + ikU + 1 + l 2 k 2 )t 22 = / 2 (1 - A)D 2 t 22 + 2ikal 2 A(D 2 - k 2 )Di\) 

+ 2k 2 aipT 12 + 2ika(T n - l)Dip (41) 

{-iuj + ikU + l + l 2 k 2 )t n = l 2 (l - A)D 2 t n - 2ikl 2 A(D 2 - k 2 )D^ 

+ (1 - a 2 )7ti2 + iktpDTu + 2ikDip + T 12 [(l - a 2 )D 2 i) - k 2 {\ + a 2 )^] (42) 

{-iuj + ikU + l + l 2 k 2 )t l2 = 1 2 {1 - A)D 2 t l2 - l 2 A(D A - k 4 )^ 

- jt n + ikijDT 12 + (1 - T n )(D 2 + k 2 )ip + 2(1 - a 2 Y l k 2 i,T lx . (43) 



5 Equivalence of the stability analyses at / — > and / = 

5. 1 Introduction 

The limit / — > is a singular limit, in the sense that the coefficient of the 
highest derivative in the governing equations is multiplied by I 2 . As discussed 
above, therefore, the one-dimensional solutions to equations (10-19) are very 
different for small / from / = 0: for I ^ the shear stress T& at which a banded 
state can exist is uniquely selected, whereas for / = banding can occur at 
any shear stress Ti < T b < T 2 . 

In this section we consider the limit of very small I in the two-dimensional 
stability problem, and show that, in contrast to the one-dimensional case, 
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the limit / — > is equivalent to / = (with suitable interfacial boundary 
conditions), for any given base state. To do this, we start by considering the 
full diffusive model, in the asymptotic limit of small I. As I — > the diffusion 
terms become unimportant everywhere in the flow apart from in the inner 
region, which has dimension /. This region separates two phases of flow at 
different shear rates, in which the governing equations are those of the original 
(non- diffusive) Johnson-Segalman equation. 

The crucial question to address here is what boundary conditions are imposed 
at the edge of each of these phases as a result of the structure of the inner 
solution. Our central finding is that these conditions are the same in the limit 
of a thin interface as they would be for the case of an infinitesimally sharp 
material interface (that is, an interface across which material does not pass) 
with I = 0. Accordingly, we now summarise the case of a sharp (diffusion-free) 
material interface / = before proceeding with our analysis of weak diffusion, 
1^0. 

McLeish [36] has argued that if the functional relating current stress to strain 
history is well-behaved, a small perturbation to the flow can only cause a small 
perturbation to the stress in any fluid element. Thus material can never be 
transported across a sharp interface since this would cause an order 1 change 
in the stress for that material. 

This physical argument leads to a kinematic boundary condition: if the in- 
terface is defined by points y = Tj + ((x,7j,t) (where Tj is the Lagrangian 
cross-channel coordinate) then 

£ = (44, 

This boundary condition is standard for an interface separating two different 
materials: in particular, Renardy [37] applied it to two phases of a shear- 
banded Johnson-Segalman fluid. For our linear perturbation problem, if the 
interface Tj = k is displaced from y — k to y = Tj + (£, then since Tj = y — (£, 
the kinematic boundary condition (correct to linear order) is expressed as 

{-iuj + ikU(K))C = -ikip(K). (45) 

Then imposing continuity of velocity and traction across the interface gives 
boundary conditions for the jumps in ip, Dip and s^- across the interface. 

With this summary of the diffusion-free sharp material interface in mind, we 
now proceed to show that the case I — > is equivalent to it. 
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5.2 Asymptotic analysis 



We now consider the asymptotic limit I — > 0. We scale our perturbation flow 
so that the values of standard quantities {e.g. ip, Dip, t 12 ) are order 1 in the 
outer regions. Within the inner region we pose the following asymptotic series 
(in which superscripts on ip, s^, tij and p are for labelling only, while those 
on I indicate powers): 

i> ~ ip° + lip 1 + i 2 ip 2 + • • • Uj ~ r 1 *-. 1 + 4 + • • • (46) 
p ~ rV 1 + p° + ■ ■ ■ ~ r 1 *-. 1 + 4 + • • • (47) 

These scalings are not obvious a priori; rather, the existence of a solution 
with these scalings will justify its choice. As before, we scale lengths within 
the interfacial region as £ = (y — k)/1, giving ID = d. This analysis will be 
valid as long as this lengthscale is well separated from other lengthscales in 
the flow: in particular, the arguments and analysis below will not apply for 
very short waves for which k ~ I -1 : we restrict our analysis to the case hi <C 1. 

We also pose a series for the eigenvalue u: 

uj ~ ujq + lio\ + • • • (48) 

Well away from the interface, all quantities are at most order 1, which we have 
ensured by scaling the whole perturbation. This gives a series of conditions on 
our variables as £ — > ±oo. In particular, since Dip, Uj and s^- must be finite, 
with no contribution at 0{l~ r ), we have the conditions 

dip -> 0, fr. 1 -> and s"- 1 -> as £ -> ±oo. (49) 

We split the governing equations into successive orders of I. At order l~ 2 we 
have, from the stress equations, simply (from (39)) 

d 2 ip° = (50) 

which, along with the matching condition (49) gives ip° = constant = cto- This 
means that the streamfunction ip (or velocity in the y-direction) is continuous 
across the interface region, and a is the value of ip at the interface. 

At 0(l~ 2 ) in the momentum equations, we have: 

ds^ = ds^i = (51) 

which, along with (49), give = s 2 i = 0. Then the stress equations (38)- 
(40) at O^- 1 ) give 

p- 1 = -(1 + a)-H^, = 2(1 - a 2 )-H£, (52) 



13 



edV + tii = 0. (53) 
The equations governing t" 1 at order 0(/ _1 ) become: 

{-iujQ + ikU{n) + l)t n 1 = (1 - A)^" 1 + (1 - a 2 )^ 1 

+ ika dT u + (1 - a 2 )T l2 d 2 ^ 1 (54) 

(-iw + ifcl7(K) + l)^ 1 = (1 - A)rf 2 tr 2 1 - AdV - 

+ ika dT 12 + (1 - 7n)dV (55) 

* = 0. (56) 

We now introduce a quantity 

£ = — ikao/(— iouo + ikU(K,)). (57) 
Equations (53)-(55) then give the following system: 

ee*V + tii = (58) 

t£ = (1 - A)rf 2 t n 1 + (1 - a 2 )ytn + (1 - a 2 )T 12 d 2 ^ 

- (-iw + ifcl/fa))^ 1 + C^n] (59) 

t^ 1 = (1 - A)rf 2 t 12 1 - AdV 1 - 7^n + (1 - 7n)dV 

- (-iw + ifcE^/c))^ 1 + C^T 12 ] (60) 

which is solved by 

dV = -Cd-y, t'12 = ~CdT 12 , t£ = -(dT n , (61) 

all of which tend to zero for large |£| as required. These are the first deriva- 
tives (in the flow-gradient direction) of the corresponding base-state quanti- 
ties, which for long waves are the perturbations we would expect from a simple 
displacement of the interface by an amount (. Integrating d 2 i\) x gives 

d^ 1 = e*i - C7- (62) 

We note that (52) also gives 

= 2(1 - a 2 y l t^ = -2(1 - a 2 )- 1 ^!! = ~C^n- (63) 
At 0(/ _1 ), the two momentum equations give 

iks^ + ds° l2 = iks^ + ds° 22 = (64) 
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and thus 



s? 2 = oi2 + ik(Sn 



S°22 = "3- 



(65) 



5.3 Physical meaning 

We have now calculated all the singular components of the perturbation flow 
for the asymptotic limit I — > 0. What conclusions can we draw about the 
behaviour of the outer quantities (which are not directly affected by diffusion) 
at the edges of the inner layer? Looking at velocities and stresses (and recalling 
that Dip ~ l~ 1 dip so that the order 1 term in Dip comes from ip 1 ), we have 
the following conditions (to within corrections of 0(1)): 

i/)(k-) = i/>(k+) = a (66) 

which follows from (50), 

+ CV(«-) = + Ct(«+) = «i (67) 

from (61), and two conditions from (65): 

Si2(«-) - ik(Sn(K^) = si 2 (k+) - ik(S n (K, + ) = a 2 (68) 

S22(K-) = S22(«+) = 013 (69) 

along with the definition 

( = — ikao/(— iujq + ikU(n)) (—iuo + ikU(K,))( = —ikip(K,). (70) 

We have used the notation X(k-) to denote the value of variable X in the 
limit y — ► k of the outer region given by < y < k; similarly, X(k+) denotes 
the value as y — > k in the outer region k < y < 1. 

These equations exactly reproduce those for a sharp material interface with 
displacement ( from its original position and a discontinuity in both 7 and 
S11 across it in the base state. If we identify the quantity ( as a displacement 
in the interface, then the conditions (66-67) correspond to continuity of the 
total fluid velocity at the interface, and (68-69) come from the momentum 
equations, and correspond to continuity of the total traction (g + g_£) ■ n at 
the interface, where n is the normal to the displaced interface. Of course the 
interface itself has width of order /, so only in the asymptotic limit I — > 
do we truly have an interface to consider, which is why the analysis above is 
only valid in this limit. Note that the original Johnson- Segalman equations 
governing a linear perturbation are fourth-order in ip, and equations (66-67) 
supply four boundary conditions at the interface, as we require. 
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In this limit I = but with a finite interface displacement, we have effectively 
exceeded the validity of our original linearisation. When we linearised about 
the base state, we implicitly assumed that the scale of all perturbation ve- 
locities, displacements and so on was much smaller than any other quantity 
in the problem, and in particular that displacements must be smaller than 
the width of the interface region, I. In taking the asymptotic limit / — > we 
have continued the assumption that displacements are small relative to the 
interface size: so this asymptotic limit is not founded on a well-posed physical 
problem. This makes it even more surprising that the limit should correspond 
exactly to the physical problem with an infinitely sharp material interface, in 
which the (small) displacement is much larger than the (zero) interface width. 

If we consider a perturbation to consist of an initial perturbation of the inter- 
face (this assumption is not always valid, as demonstrated recently by Kupfer- 
man [40]), then the driving terms for the perturbation are the ( terms in (67) 
and (68), which depend on discontinuities in the base-state. If the shear rate 
were continuous across the interface, the only driving term would be the dis- 
continuity in Su in (68); this would then allow us to recapture the normal 
stress instability mechanism elucidated by Hinch et al. [2]. The discontinuity 
in shear rate in (67) provides another, more complex, mechanism for instabil- 
ity, which is not yet fully understood. 

At order /, however, equations (66-69) will not be satisfied so the diffusion does 
alter the zero-transport property. This is to be expected if the stress diffusion 
derives physically from diffusive transport of individual polymer molecules 
carrying their stress history with them: we would expect material to travel 
across the interface, a distance of /, with diffusive velocity of I 2 , on a timescale 
of order I -1 . This timescale corresponds to a contribution of order / to cu, i.e. 
the next term in the asymptotic series, lo\. 



5.4 Wall boundary conditions 

In our discussion of the boundary conditions applying across the diffusive 
interface layer, we have implicitly assumed that this layer is the only place 
in the perturbation flow where diffusion is important. However, the stress 
boundary conditions Do~ij = at the walls are not automatically satisfied by 
the perturbation flow of a pure Johnson-Segalman fluid. We would expect that 
a stress boundary layer would form at each wall, matching between our outer 
solution and the no-flux boundary conditions. 

This does indeed occur, but the influence of these stress boundary layer regions 
does not extend into the bulk of the flow, and they do not affect the stability 
of the flow. They are standard boundary layers as one would expect for what 
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is essentially an advection-diffusion equation in extra stress. Although we are 
not presenting an analytical justification of the assertion that their effect is 
negligible, the close match between our results (neglecting the wall boundary 
layers) and those of [1] and of §7.2 at small finite I with no such assumptions, 
justifies our claim. 



6 Long wave analysis <^C 1 at I = 

We showed in Section 5 that as long as the wavenumber k is not very large 
(we assumed in Section 5.2 that kl <C 1), the interface between the two phases 
of the flow may be treated as a material interface, with corrections appearing 
at 0(1) for small /. A typical estimate of / (the ratio of mesh size to channel 
width in [41]) is 2 x 10~ 4 , so the waves which are excluded from our analysis 
are very short indeed. 

In summary, for small diffusion, the diffusive terms in the outer region do 
not affect the flow at any physically important level, so the only role played 
by diffusion is in the "boundary" conditions that it imposes on the outer at 
the edges of the inner, interface, region. We have just shown that these are 
equivalent to treating the interface as a material interface in the case I = 0: 
so the solution to the I = problem with a material interface must be the 
same as the I — * limit of the diffusive problem. This equivalence is verified 
in Section 7 by the agreement between the / = results of this paper with 
earlier numerical results by one of the authors [1] for small /. 

This finding allows us now to consider the simpler (although artificial) limit 
I — 0, in which the diffusive terms play no part at all, and the base state 
quantities 7, T n and T 12 are constant within each fluid phase, and treat the 
interface as a material surface. We are now dealing with the original Johnson- 
Segalman fluid [19] and considering a stability problem which was studied 
by Renardy [37] in 1995, although we have some analytical results for long 
waves where her work was purely numerical. However, as discussed above, for 
planar Couette flow of a Johnson-Segalman fluid in the non-monotonic region 
of the flow curve, there is no mechanism for selection of the shear stress. 
Any shear stress that intersects both the low- and high-shear-rate parts of 
the flow curve, T& with T\ < Tb < T2 (figure 1) is a possible solution to the 
governing equations, with the choice of shear stress (for a given wall velocity) 
governing the position of the interface. Renardy's solution to this difficulty 
was to assume top-jumping, that is, the maximum value of the shear stress for 
which the flow curve is multi-valued. We, however, are using this local model 
as an approximation to the non-local case / 7^ 0, so we adopt the unique value 
of Tb selected by use of the gradient terms in the stress-diffusive non-local 
case. 
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In the case of long waves, k — > 0, p, Sn and S22 are larger than the other 
perturbation quantities by a factor of fc -1 , and the eigenvalue u is expected to 
be of order k, based on prior experience of long- wave interfacial instabilities [2]. 
If we scale these quantities accordingly: p = k~ l p, Sa = fc -1 Sjj, uj = ku, then 
the governing equations within each fluid, from (37-43), become, for the force 



balance: 

isu + Ds 12 = ik 2 s 12 + Ds 22 = 0, (71) 

for the total stress: 

sn = ~P + 2eik 2 D^j + (1 - a)' l k{t 22 + t n ) (72) 
S12 = e(D 2 + k 2 )ip + ti2 (73) 
S22 = -p- 2eik 2 Dij + (1 + a)- l k{t 2 2 - *n), (74) 

and for the polymer extra stress: 

(-ikZJ + ikU + 1)*22 = 2aik{T u - l)Dip + 2ak 2 ipT l2 (75) 



(-iA;57 + ifct/ + l)fn = 2ikDijj + (1 - a 2 )7t i2 

+ T 12 ((l-a 2 ) J D 2 ^-(l + a 2 )A; 2 ^) (76) 



(— ifcaJ + iA;Z7 + l)*i 2 = (1 — T n )(D 2 + k 2 )ip — jt n + 2(1 — a^^Tnk 2 ^. (77) 

If the interface has been displaced from its original position Tj = k to a new 
position fj = k + ( exp (i/cx — iut), then the conditions of continuity of velocity 
and traction across the interface give us continuity of the following quantities 
at y = k: 

-ikip iC + Di) -ik(S u + s 12 s 22 (78) 
in agreement with the results of (66-69) for I — > 0. 

Finally, the condition that the interface should be a material surface, equiva- 
lent to (70), gives 

(-u + U( K )){ = -rl>( K ). (79) 

This system of equations is essentially fourth-order in ip, so the general solution 
within one shear band will contain four unknowns. The no-slip, no-penetration 
boundary conditions on the walls remove two unknowns in each band, leaving 
a total of four unknowns to be determined through these jump conditions at 
the interface. Since there are four jump conditions here, and no driving term 
not proportional to ip, the existence of a non-zero solution to the problem is 
the eigenvalue condition which allows uj to be fixed. 

We expand each perturbation quantity as a regular power series in k. A critical 
quantity in each fluid phase is the marginal viscosity, that is, the slope of the 
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stress-shear-rate relation a(j) at our selected shear rate: 



" = "df = d^ 7 + Tl2) = 6 + (l + (l-a 2 ) 7 2 ) - (80) 

Using to denote the value of quantity X in the lower shear-rate phase, and 
X H to denote its value in the higher shear-rate phase, the solution at leading 
order in k is 

S12 = JL L D 2 ^ L = J[ H D 2 ^ H = Ay + B s u = s 22 = -p = iA (81) 

where A and B are constants. The physical interpretation of this result is that, 
to this order, the perturbation consists of simply displacing the interface: a 
one-dimensional perturbation. If uo = (as it is to this order in k), then the 
system is still on the "adiabatic" constitutive curve of figure 1, but with a 
modified stress T b . Thus the change in stress is 

si 2 = AS 12 = = /LD^S (82) 

since the change in shear rate at this order is D 2 ip. 

Integration of D 2 ip subject to zero perturbation velocity on the channel walls 
determines ip in each phase, and the interfacial continuity conditions then give 
the values of the two constants: 



A _ 6([JI h k 2 - -p L ( K - l) 2 ]]l I jI H (<y H - 7 L ) 

[Ph*- 2 - ~P l {k - l) 2 ] 2 - Ajl L -p H K(K - 1) ' 
B = 2 ([~Pl(k 3 - 3k + 2) - ]I H K, 3 ]p I jI H ('y H - 7l) ^ 
[JI h k 2 - JI l (k - l) 2 ] 2 - AJ1 l JI h k(k - 1) 

The eigenvalue at this order is real, which corresponds to the perturbation 
travelling in the flow direction, and is given by 

„ „ ^ + ^f^r™*-™ , j + Oie). (85) 

[fi H K, 2 - fi L {K - l) 2 ] 2 - 4fl L fi H K{K - 1) 

The imaginary part of u determines the stability or instability of the flow 
(Im(u;) < for stable flow, Im(cj) > for unstable flow). This appears at the 
next order. This problem can still be solved analytically, but the result is too 
unwieldy to reproduce here. 

Instead, in figure 2 we show the stability boundaries for this long-wave mode 
of perturbation, plotted in the parameter plane of T h ^J (1 — a 2 ) against k. 

It will be seen that the range of parameter values over which the flow is 
unstable includes most of the available values of k, with the exception of 
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0.2 0.4 0.6 0.8 1 

Proportion of low shear rate band, k 

Fig. 2. Stability boundaries in parameter space for long-wave perturbations (stability 
or instability at order k 2 ). The x-axis is k, the proportion of low-shear-rate fluid in 
the base flow; the y-axis is the scaled shear stress (f — a 2 ) 1 / 2 ^. The curves shown 
are for e = 0.05, and these boundaries are independent of the value of a. The solid 
horizontal lines show the minimum and maximum possible values of (1 — a 2 ) 1 / 2 !],, 
being 0.435207 and 0.552806 respectively for this value of e. The dotted horizontal 
lines show the selected values of (1 — a 2 ) 1 / 2 ^ for the two limiting diffusive models: 
Olmsted's model (A = 0) predicts (1 - a 2 ) l/2 T b = 0.48284, while Yuan's model 
(A = 1) has (f - a 2 ) l / 2 T b = 0.45 for e = 0.05. U denotes the unstable region 
(typically for two roughly equal shear bands) and S the stable regions (typically 
one or other shear band being very narrow). The stable regions here have only been 
shown to be stable to asymptotically long waves: as we see in section 7, most of 
these parameters do show instability for some value of the wavenumber, k. 

values giving an interface very close to one of the walls. In particular, for the 
shear stress which is selected by the stress-diffusion model A = 0, there is 
instability for 0.108 < k < 0.933. 

In the next section we give the results of numerical calculations at finite values 
of k, and we will show these analytical asymptotes along with those numerical 
results. 



7 Numerical Results 



In order to access perturbations with moderate or short wavelengths, we solve 
the linearised equations numerically. The scheme used is a shooting method, 
in which a value of u> is guessed and the Newton-Raphson method is used to 
find the true value. The equations are integrated inwards from each wall and 
the jump conditions at the interface provide the dispersion relation through a 
determinant condition as introduced by Ho & Denn [42]. 
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1 2 3 4 5 6 7 

Average shear rate, C/ wa n 
Fig. 3. The dependence of the growth rate of instability on the average shear rate 

U wa ii and on wavenumber, k. The parameters are e = 0.05, a = 0.3 and the shear 
stress is that selected by the non-local model with stress diffusion, T& = 0.506158. 
The bold solid curves mark the boundaries between instability and stability (the 
unstable region being the centre of the plot), the thin solid lines mark local maxima 
in the plot of growth rate against wavenumber, and the dashed curves mark local 
minima in the same plot. 

7.1 Most unstable mode 

In this section we choose as illustrative the parameters e = 0.05 and a = 
0.3, and investigate the wavenumber-dependence of the instability. We use 
the nonlocal model with stress diffusion (i.e. A = 0) to select the shear 
stress T b = 0.506158 of the one-dimensional base state, but the subsequent 
two-dimensional linear stability calculations are carried out using the local 
Johnson-Segalman model. 

We vary the average shear rate across the channel, £/ wa n, which changes the 
proportion of the lower shear rate band. Our nondimensionalisation means 
that C/ wa u is the same as the Weissenberg number. The behaviour of the in- 
stability against wavenumber of the perturbation changes as C/ wa u varies. It 
is only for extremely low shear-rates, that is, thin high-shear bands that this 
flow is stable for all wavenumbers. In figure 3 we plot various curves in the 
U wa \\-k plane. The thick solid curves mark the borderline between unstable 
and stable parameters. For most values of C/ wa n in the shear-banding region 
0.661 < C/ wa n < 7.089 these curves do not appear as the perturbations at all 
wavenumbers are unstable. 

The thin solid curves gives the wavenumbers at which the growth rate reaches 
a local peak. Where there is only one such curve (for example, at C/ wa u = 4), 
it is tempting to regard this as the most unstable mode; however, within the 
confines of our local approximation to the diffusive Johnson-Segalman fluid, 
there is an instability to very short waves [37], with a growth rate independent 
of C/ wa n, and for many layer arrangements this is in fact the most unstable 
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mode. For our illustrative parameters, the growth rate of the k — > oo mode is 
0.533, and this mode is the most unstable for U wa \\ > 3.2. For lower values of 
U wa n the most unstable mode is the peak at the lowest wavenumber, although 
this wavenumber also becomes very large as C/ wa n — > 7^. In practice, as we 
shall see in §7.2, for realistic non-zero values of I we expect diffusive effects 
to damp very short-wave instabilities, and the peak growth rate will always 
occur at finite k: in §7.2.1 we will show how the maximum growth rate varies 
with average shear rate, with diffusion added in a semi-empirical way. 

In figure 4 we give four plots of growth rate against wavenumber, to illustrate 
the different types of behaviour seen in figure 3. The parameters chosen are 
U wa \i = 1, for which long waves are stable, and instability arises at finite 
k] U wa u = 2, which has one maximum which is the most unstable mode; 
U wa \\ = 4, which has a single local maximum but for which the most unstable 
mode is for asymptotically short waves; and U wa \\ = 6.8, for which (again) 
long waves are stable, and two separate unstable regions 0.51 < k < 2.1 and 
k > 5.45 are observed. In each case the behaviour as k — > 00 is the short-wave 
instability predicted by Renardy [37], with growth rate 0.533, and the long- 
wave behaviour matches the analytic calculation of section 6. The long- wave 
asymptotes are plotted along with the numerical calculations. 



7.2 Comparison with results at small finite I; Fielding (2005) 

In this section we carry out a comparison with the numerical study of the 
full diffusive Johnson-Segalman model at A = 0, first published in [1]. From 
the findings of section 5.2, we expect that, for a given one-dimensional base 
state, the numerical results of [1] should converge, in the limit I — > 0, to the 
asymptotic results calculated here at I — 0. 

In figure 5, therefore, we reproduce (as points) the data from Figure 3 of [1] 
showing instability for small finite values of I, along with (as curves) numerical 
results for the pure, I = Johnson-Segalman fluid under the assumption that 
the interface is a material surface, and the long-wave asymptotic form for 
the same / = situation. On the left we present the raw data. From the 
numerical data we observe that Im(cj) ~ ak 2 + b(l) for long waves k < 1, with 
b(l) PS —10/. We confirm the scaling of b(l) with I in the appendix. On the 
right, therefore, we add an additional term of 10/ to the numerical results for 
finite /, and observe that for long waves this collapses all the points onto the 
/ = curve. 

The numerical factor 10 which is used at order / for the long-wave results is 
not calculated analytically, but deduced from the k — > intercepts from the 
data for the different values of /. However, we can make some progress towards 
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Fig. 4. Plots of growth rate against wavenumber for four different values of £/ wa u, 
the average shear rate. The parameters here are e = 0.05, a = 0.3, I = 0, and 
TJ, = 0.506158, the value predicted by the stress-diffusion model. In each case the 
order k 2 behaviour of long waves is a dotted curve, (a) £/ wa ii = 1. Long waves 
are stable with decay rate Im(u;) ~ — 0.0060/c 2 and instability begins at k > 0.65 
(not discernible on the scale of the plot). The region k > 30, in which the growth 
rate tends smoothly to 0.533 from above, is omitted to allow a clearer view of the 
behaviour for longer waves, (b) U w& \\ = 2. Long waves are unstable with growth rate 
Im(u;) ~ 0.15/c 2 , and the most unstable mode is at k = 4.3 with growth rate 0.698. 
(c) ?7 wa ii = 4. Long waves are unstable with growth rate Im(u>) ~ 0.42/c 2 , and there 
is a peak in growth rate at k = 2.3, but the most unstable mode is k — > oo. (d) 
U wa ,n = 6.8. Here, as in (a), asymptotically long waves are stable. Long waves have 
decay rate Im(u;) ~ — 0.00076/c 2 . We have omitted the region k > 30, in which the 
growth rate tends smoothly to 0.533 from below. 

calculating this value. A full derivation of this process is given in appendix A: 
the final conclusion is that, as k — > and I — > 0, the growth rate can be 
written asw = ia l + 0{l 2 ) +0(k), where 



0o 



4[h h k 3 - (i L ( K - l) 3 ]/i L /i^(7^ - 7 L ) 



a 



[fi H K 2 - fi L ( K - l) 2 ] 2 - 4{i l {i h k(k - 1) 



(86) 



and 5" is an unknown parameter which depends on a, e and T& but not on 
C/ wa n or k. For the values of a, e and T b associated with figure 5 the data 
suggest that a ~ 14. In figure 6 we plot a numerical calculation of cr against 
k for these values of a, e and T b . The curve is given by the prediction of (86) 
with a = 14. The agreement between theory and numerical calculation of this 
term is remarkable. As k — > 1 our prediction (with a = 14) is a pa —130: 
numerically it is not possible to investigate extremely narrow bands as the 
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Fig. 5. Growth rate Ym.(u) plotted against wavenumber, k for Couette flow at 
e = 0.05, a = 0.3, ?7 wa ii = 2, Tb = 0.506158. Curves: / = 0, calculation 
for a non-diffusive Johnson-Segalman fluid with no material transport across the 
interface between phases. The dashed curve is the long-wave asymptotic form 
Im(u;) ~ 0.1506/c 2 ; the solid curve is the full numerical calculation. Points: (from 
highest to lowest) / = 0.00125, 0.0025, 0.005, 0.01. The figure on the left shows the 
true growth rate; on the right for the finite-/ results we have plotted lm(uj) + 10L 




U 0.2 0.4 0.6 0.8 1 

Proportion of low shear rate band, k 
Fig. 6. Dependence of the coefficient of I in the growth rate of very long- wave 
perturbations on the interface position, k. Here a = 0.3, e = 0.05 and Tb = 0.506158. 
The points are from numerical calculations and the curve is given by the analytical 
prediction from (86), with a chosen to match the data at k = 0.79. 

diffusive layer needs to be clear of the walls, but the numerical calculations of 
this long wave intercept match the theoretical prediction well even at k — 0.97 
where a « —48. 

This diffusion-induced stability to very long waves k — > was to be expected. 
Calculations had already shown that an interface between shear bands at this 
selected shear stress should be stable to one-dimensional perturbations: that 
is, if the whole interface is rigidly displaced from the selected position it should 
relax back there. A simple displacement of the interface corresponds to the 
limit of very long waves, and so this exponential relaxation is precisely the 
negative growth rate we have calculated in the limit k — ■> with / small but 
finite. 



24 



0.6 




+ + + + 



+ 



x 



+ 



X 



+ 



X 



+ + 



X 



X 



+ 



X - 



+ 



-0.2 - 



+ 



0.3 



i i i i J: i i i i i i 1 i i 

1 23456789 1011 12131415 

Wavenumber, k 



Fig. 7. Dependence of the growth rate of the instability on wavenumber. Parameters 
e = 0.05, a = 0.3, T& = 0.506158, f/ wa ii = 4. The long-wave mode leads to a peak 
growth rate at k « 2, but the growth rate increases again for shorter waves, and 
the most unstable mode (in the limit I = 0) is for very short waves. The solid curve 
is for / = 0; the points (from highest to lowest) are / = 0.00125, 0.0025, 0.005, 0.01. 

In figure 7 we show another comparison between the I = and / ^ 0, A = 
cases, this time at C/ wa u = 4. The convergence of the / ^ results to the 
I = case as I — > is clearly visible. At I = the growth rate for long waves 
is Im(cj) ~ 0.42/c 2 , and the calculation of section 7.2 predicts a = —10.9 
for these parameters, so for long waves at small finite / we expect Im(u;) ~ 
-10.9/ + 0.42A; 2 , which means that instability first appears at k « 5.1/ 1/2 . This 
scaling k ~ l 1 ^ 2 for the lowest unstable wavenumber is universal provided that 
long waves are unstable in the I = case. 

For very short waves, we can see from figure 7 that the addition of diffusion 
has a large effect on the eigenvalue. We expected this when we stated in 
section 5 that our analysis would only be valid for k <C I -1 . In fact we can 
see empirically that for these parameters, Im(cj) ~ —Ylkl + 0(1) for fixed / 
as k — > oo, leading us to predict instability for k / _1 and stability for very 
large k. The size of the pref actor in this case means that the results for finite 
/ deviate from the I = limit earlier (as k increases) than a simple scaling 
argument might have led us to expect. 

In summary, if 0.108 < k < 0.933 (for the stress-diffusion model) then in 
the limit I = the flow is unstable to perturbations of all wavenumbers, and 
in this case we expect the diffusive flow to be unstable over a large region 



As a final comparison between our calculations and the numerics for small 
finite /, in figure 8 we give the perturbation streamfunction ip and its first 
derivative (proportional respectively to the cross- and along-channel velocity 
components) for one specific mode. We have plotted the real part of ip and its 
derivative at the phase (x-position) where i/)(k) is real. The parameters (given 



/ 1/2 <A;</- 1 . 
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Channel position, y 
Fig. 8. The streamfunction ip, proportional to v (continuous function) and its deriva- 
tive Dip, or u (discontinuous) of the least stable mode at e = 0.05, a = 0.3, J7 wa ii = 2, 
T b = 0.506158, k = 0.1. This mode is unstable at I = 0, with growth rate 0.00148; 
at I = 0.005 it is stable with growth rate -0.0485454. The crosses are from the full 
calculation at / = 0.005, and the solid lines (lying under the crosses except near the 
interface) from the limit I = using the material surface condition. The stream- 
function is plotted against position across the channel, y, and is normalised such 
that iP(k) = 1. We show only the real part of ip and Dip here; the imaginary part is 
smaller by a factor of order k. 

in the caption to figure 8) are such that, in the limit I = 0, long waves are 
unstable and the mode k — 0.1 shows this growth. However, for / = 0.005 the 
decay term at 0(1) dominates and the mode is stable. Nonetheless, the form 
of the streamfunction (equivalent to lOi times the y-velocity in this case) and 
its first derivative (which is the x-velocity) is extremely similar between the 
two modes. 



7.2.1 Most unstable mode with diffusion 

Finally, in figure 9 we give an empirical idea of the maximum growth rates 
which might be expected to be seen for a range of average shear rates, and for 
two different values of the slip parameter a. For the pure JS model, the most 
unstable mode is often the short-wave limit k — > oo, which will be stabilised 
by diffusion for the modified model. It is therefore unhelpful to give the growth 
rate of the most unstable mode; rather, we attempt to give the expected most 
unstable growth rate in the presence of some diffusion I = 0.00125. In order 
to carry out the computations for a large variety of parameters, we make two 
very broad assumptions: 

• Based on curve fitting of the data in figure 7, the growth rate at any k and 
any a may be reasonably approximated by 

Im(^) approx = Im(u;)js + o~ l - 13kl. 
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Fig. 9. Maximum growth rate plotted against average shear rate for a semi-empirical 
model of the stress-diffusion fluid with I = 0.00125, e = 0.05, and A = 0. Lower 
curve a = 0.3; upper curve a = 0.8. In both cases regimes with very narrow bands 
of one or other shear rate are stable to all perturbations. The wavenumber of the 
most unstable waves increases as we approach these stable regions. 

• In determining <t , the value a in equation (86), which is a function of a, 
e and Tf,, will be taken to be 14 (the true value at a = 0.3, e = 0.05 and 
T b = 0.505158) independent of these parameters. 

Using these two assumptions and the numerical calculations for the pure JS 
problem, we predict the wavenumber and growth rate of the most unstable 
mode for two different fluids. We use the stress-diffusion model to select the 
shear stress in each case, and consider the cases e = 0.05, a = 0.3 and e = 0.05, 
a = 0.8. 

There are configurations which are stable to all perturbations, which are 
those layer arrangements with one or other shear band being very narrow, 
i.e. U wa \i ~ 7l or LV" wall rj j H , However, the vast majority of mean shear-rates 
in the shear banding regime produce banded flows which are linearly unstable 
with a moderate growth rate. The figure suggests that influence of the larger 
slip parameter a = 0.8 tends to enhance the instability, but we cannot draw 
concrete conclusions from such an empirical model. 



7. 3 Comparison with other previous studies 
7.3.1 McLeish (1987) 

An early paper by McLeish [36] considered capillary flow with a constitutive 
equation with a non-monotonic flow curve. He predicted exactly the opposite 
of the long- wave behaviour we have found: for slow {i.e. long- wave) pertur- 
bations, he found instability only for very narrow regions of high-shear-rate 
material (the band close to the wall). That work used a slightly different con- 
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stitutive equation [43], based on reptation theory for linear polymers, but his 
model has similar behaviour to ours in steady shear. He described the stability 
property of the flow as a dependence of the throughput in the lower-shear rate 
region on both the absolute position and the slope of the interface, with a for- 
mulation which is mathematically equivalent to ours except for the different 
base flow. 

In terms of instability mechanism, McLeish suggests that since a normal stress 
difference is required in order for the system to "see" the gradient of the in- 
terface perturbation, normal stress effects are critical to the instability. From 
our calculations we see that this does appear to be true, but that the mech- 
anism of instability is not quite the clean recirculation mechanism found by 
Hinch et al. [2] for coextruded fluids having matched viscosities and a jump in 
Aq across the interface between them. In our equations there are two driving 
terms: the jump in Dip, proportional to the difference in the base-state shear 
rate across the interface and the interface displacement (67); and the jump in 
S12, proportional to the difference in the base-state Aq across the interface and 
the slope of the interface (68). Algebraically, we can artificially separate these 
out, and in most cases studied here the normal stress term was weakly sta- 
bilising, and the instability comes from the interaction of the shear-rate-jump 
term with the normal stresses in the bulk of each fluid. 



7.3.2 Renardy (1995) 

Renardy [37] examined the stability of the local JS model in planar banded 
Couette flow. She found linear instability for short wavelengths (wavenum- 
ber greater than 8). For mainly historical reasons, however, she happened to 
confine her study to a base state corresponding to "top-jumping" (Tf, = T 2 ). 

As a check on both our analysis and our numerical eigenvalue calculation, 
we reproduce figure 2 that paper. In Figure 10 we show our own numerical 
calculation, which duplicates her results as far as can be seen from the graph 
in [37]. In the second part of the figure we plot the same growth rate again, 
along with the long-wave asymptotic form for the growth rate as it depends 
on wavenumber. We have restricted the scale in this second graph in order to 
better see the accuracy of the long- wave result. 

Since neither of the two non-local models predicts the top-jumping stress T 2 
as the selected stress, this work is unlikely to be directly relevant physically; 
moreover, the bulk of [37] focuses on a short-wave instability. While this insta- 
bility does occur for the original JS model systems at the true stress selected 
by a nonlocal model (for instance, at e = 0.05 and a = 0.3 using Tj, as se- 
lected by stress diffusion, it has growth rate 0.533), when diffusion is added 
the short-wave instability mechanism is destroyed and this mode will not be 
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Fig. 10. Reproduction of figure 2 from Renardy [37]. Parameters are e = 0.05, 
a = 0.8, Tb = 0.921. Growth rate is plotted against wavenumber. The upper fig- 
ure reproduces the original exactly (to the naked eye; the original data were not 
available) and in the lower figure we reproduce the long-wave portion of the graph. 
Solid curve: numerical calculations; dotted line: asymptote — 0.0312/c 2 from long 
wave calculation. As stated in the text of [37], instability appears for waves shorter 
than k 8, and the peak growth rate of instability is for waves having k ~ 41. 

seen. 



7.3.3 Yuan (1999) 

In 1999, Yuan [39] carried out time- dependent simulations using a model cor- 
responding to A = 1 in our model. He found that, for e = 0.05 and a = 0.8, 
the system uniquely selected a shear stress of T& = 0.81 ±0.04 for any aver- 
age shear rate in the range over which a homogeneous solution would be on 
the unphysical descending branch of the constitutive curve. This gives values 
T 6 * = (1 - a 2 fl 2 T b = 0.486 ± 0.024, which is slightly outside the true value of 
0.45 from equation (31). 

He found stable steady states for all parameters in this range. Our long- 
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wave analysis predicts instability to long waves at his parameter values for 
0.186 < k < 0.974, that is, 1.13 < U wa \\ < 7.56, suggesting that several of his 
simulated flows (at C/ wa n = 2, 3, 4, 5 and 6) should be unstable to long- wave 
perturbations. However, he was simulating a finite length of channel (and with 
a relatively large value of /) which allows access only to specific wavenumbers, 
which may be the reason that he did not capture the instability we have 
demonstrated. Alternatively, Yuan's results may be an artefact of averaging 
over time t and/or flow location, x. Note that if these simulations are time- 
averages of a nonlinearly fluctuating interface where a flat interface is linearly 
unstable, then we would expect the average shear stress reported to be higher 
than the selected value for the ID base state, as observed above. 

Another comment which seems unusual is that Yuan states that even when 
1 = 0, the interface has finite width. Although there are physical reasons why 
this may be true in practice, it is not predicted by the governing equations in 
the limit / = 0: perhaps his observation is a grid-scale effect. 



8 Conclusions and Discussion 

We have investigated the two-dimensional linear stability of plane Couette flow 
of a shear-banding fluid. If the equations governing the Johnson-Segalman 
model are regularised using a small amount of stress diffusion to provide a 
uniquely selected one-dimensional banded base state, we have shown that as 
far as two-dimensional linear stability is concerned, the limit of no diffusion is 
a regular limit in which the interface region becomes a strict material interface. 

Using this limit, we have demonstrated that there is a long wave instability 
for almost all possible positions of a shear band. Only the case of a very nar- 
row "spurted" region of high shear rate is stable to long-wave perturbations. 
This is in contrast to earlier work by McLeish [36], who found (for Poiseuille 
flow and slightly different constitutive assumptions) stability except for the 
case of a narrow high-shear rate band; and simulations by Yuan [39], which 
predict a steady interface between two shear bands in this situation. However, 
we agree quantitatively with results of Renardy [37] , who happened to look at 
a narrow region of high shear and found stability to long waves (although the 
paper focuses on a short-wave instability whose mechanism is likely to be af- 
fected by diffusion) . Our results are in full agreement with numerical stability 
calculations including diffusion terms [1]. For typical physical parameters, for 
which there is instability to perturbations of all wavelengths in the absence of 
diffusion, we have identified the scalings at which diffusion affects the instabil- 
ity. Small diffusion on a dimensionless lengthscale / will restabilise very long 
and very short waves, leaving the flow unstable to perturbations of moderate 
dimensionless wavenumber Z 1 / 2 < k <C 



30 



We have identified two driving forces for the instability: discontinuity of shear 
rate and of normal stress across an interface. The interplay between these 
mechanisms, even for long waves, is not fully understood. 

This widespread instability suggests that the existing theoretical picture of 
two stable shear bands separated by a steady interface needs further thought. 
Indeed, this result is consistent with accumulating evidence for erratic fluctu- 
ations [11,9,12,44,45] in several different shear banding systems. 

Future work will investigate the behaviour of the interface in the nonlinear 
regime, beyond the validity of this linear study. One possibility is that the 
instability saturates at a small but finite amplitude — indeed, our preliminary 
investigations suggest that this is the case. This would be consistent with a 
narrowly localised but still unsteady interface, which might be interpreted as 
steady in experiments that did not have high spatial resolutions. This might 
even reconcile early data showing apparently steady interfaces with recent 
work revealing fluctuations. 

However, if this is not the case, then the use of a Johnson-Segalman type 
constitutive model, with or without stress diffusion, can almost never produce 
agreement with any steady banded structure observed experimentally. One 
would then need a new theoretical picture to incorporate the observed shear- 
banding effects within a stable flow which is either steady or undergoes only 
small-amplitude oscillations. 
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A Calculation of the order I contribution to u in the long-wave 
limit 

In this section we derive a scaling form for the (negative) growth rate of a 
perturbation at k — for small /, in the case A = 0, i.e. the case in which the 
stress diffusion is added through diffusion of the polymer extra stress term. 
This negative growth appears at order I. 

Let us return to the full governing equations (37)-(43). We will first scale with 
k and then with /. 

Using the long- wave scalings of section 6, but allowing uj to remain order 1 
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(with respect to k), and neglecting terms of order k yields the system 

isu + Ds 12 = Ds 22 = (A.l) 



sn = s 2 2 = -p S12 = eD 2 ip + tu (A. 2) 

{-iuj + l)f n = l 2 D 2 t u + (1 - a 2 ) 7 t 12 + (1 - a 2 )T 12 D 2 ^ (A.3) 

(-iw + l)t 12 = l 2 D 2 t 12 - ytu + (1 " T n )D 2 if> (A.4) 
As before, we solve the momentum equations to have 

S12 = + -B sn = s 2 2 — —p — iA. (A. 5) 

We also scale the eigenvalue: 00 = ila + 0(l 2 ). 

We now divide the flow into three regions: the two "outer" regions where the 
base state shear rate and stresses are constant, and the "inner" region where 
base state quantities have derivatives of order I -1 . We denote by X L the value 
of a quantity in the low-shear-rate band near the wall at y — 0, and by Xh 
its value in the high-shear-rate band near the wall at y — 1. As in section 6, 
we will use the marginal viscosity 

In each outer region, where derivatives are order 1, we neglect terms of order 
/ and solve to have 

4.L _ [Ay + B](JI L - e) ^ L _ 2[Ay + B}(1 - a 2 )T^ 2 

12 ~~ — f ll — — N 1 ( 1 2V21 l^-'J 

L» 2 ^l = [Ay + S]/7Z L = [Ay 3 + 3Sy 2 ]/6/J L (A.8) 

in the low-shear band, and in the high-shear band, 

H = [Ay + B](jc H - e) H = 2[Ay + B](l - a 2 )T« 

JI H 11 Jl H [l + (1 - a 2 ) 7 |] 1 ' ) 

D 2 ^ H = [Ay + B]/JI H i, H = [A(y 3 -3y + 2) + 3B(y - 1) 2 }/QJI H . (A.10) 

Within the inner region, we scale lengths as £ = (y — k)/1 so that ID = d and 
pose the series 

t ij ^i-H^ 1 + t° ij + --- + + + (A.ll) 
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Matching these quantities to the outer solutions for large |£| yields the follow- 
ing conditions: 

t£ -> 0, dip -> 0, (Pip 1 -> as f -> ±oo; (A. 12) 



as £ — > — oo; and as £ — > oo, 

t-j - *£(«), ^° - ^(«), - £>^h(«), dV - £>Vh(«)- (A.14) 

We collect orders of / in the resulting equations. At order Z~ 2 we have simply 

dV° = (A.15) 
which, along with (A. 12) gives ip° = a . At order l^ 1 we have 

edV + tr2 1 = (A. 16) 

d^n 1 " *n + (! - « 2 )7^2 + (1 - a 2 )T 12 d 2 ^ = (A. 17) 

d 2 *r2 1 - t\2 ~ 7*n + (! - Tu)d 2 i) 1 = (A. 18) 

which are solved by 

t" 1 = -<JdTn t^ 1 = -5rfTi 2 dV = -Sd'y. (A. 19) 

Substituting the matching conditions from (A.13)-(A.14) and the outer solu- 
tions leads to the condition 

Ak + b= * s \Ph k3 ~ ~Pl( k ~ ^f^Lpniju - 7l) , a 2Q , 

[p H K, 2 -JI l (k-1) 2 ] 2 -AJI l JI h k(k-1)' 

Finally, at order 1 our remaining equations are 

ed 2 i) 2 + t° 12 = Ak + B (A.21) 

d 2 t° n - t° n + (1 - a 2 ) 7 t° 2 + (1 - a 2 )T 12 d 2 ip 2 = -a 5dT u (A.22) 

d 2 t° 12 - t° 12 - 7t?x + (1 - 7n)dV = -a 5dT 12 (A.23) 

Where equations (A.16)-(A.18) were homogeneous, these are inhomogeneous 
ODEs with a forcing on the RHS which comes from the previous order cal- 
culation. Together with the conditions on from (A.13) and (A.14), they 
provide a constraint on a . 
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Although we cannot solve this problem, we can reduce its complexity. If we 
define 

a = d 2 ij 2 /[An + B] Uj = ^/[Ak + B] a = -a 5/[AK + B] (A.24) 
then the governing equations become 

ea + t 12 = 1 (A.25) 

d 2 t n - t u + (1 - a 2 )^t l2 + (1 - a 2 )T l2 a = adT u (A.26) 

d 2 i 12 - t 12 - Tin + (1 - Tn)a = adT 12 (A.27) 
and the matching conditions, 

2(1 - a 2 )T,\ ~ Jlr - e /k . 

tn - _ h \ n ' I - 21 t 12 - ^t— as ^ -oo A.28 

2(1 - a 2 )T^ ~ /7 H - e 

* n ^ - n -L (i -21 *i2 - ^ as oo. A.29 

/%[! + (! -a 2 )7H] A* 
Neither the governing equations nor the matching conditions have any depen- 
dence on k, the position of the interface between the two shear rate bands. 
Thus there is a single value of a for each set of parameters {a, e, Tf,} indepen- 
dent of C/ wa n and k. The long-wave term of the eigenvalue uj is then given by 



-4i/[/i g /t 3 -/i L (/t-l) 3 ]/i L /i g ( 7H -7 L )a 2 
" bV-^-VF-^^-l) +0(J)+0(*). (A.30) 



For the results given in figure 3 of [1], we had a = 0.3, e = 0.05 and TJ, = 
0.506158, giving shear rates of 7l = 0.66143 and •jn = 7.0893. The interface 
position was 0.79176 and a fit of the data gave uj ~ —10/. We can deduce that 

a(a = 0.03, e = 0.05, T b = 0.506158) w 14. (A.31) 
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